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Abstract. Systems of interacting random replicators are studied using generating 
functional techniques. While replica analyses of such models are limited to systems 
with symmetric couplings, dynamical approaches as presented here allow specifically 
to address cases with asymmetric interactions where there is no Lyapunov function 
governing the dynamics. We here focus on replicator models with Gaussian couplings 
of general symmetry between p > 2 species, and discuss how an effective description 
of the dynamics can be derived in terms of a single-species process. Upon making a 
fixed point ansatz persistent order parameters in the ergodic stationary states can be 
extracted from this process, and different types of phase transitions can be identified 
and related to each other. We discuss the effects of asymmetry in the couplings on 
the order parameters and the phase behaviour for p = 2 and p = 3. Numerical 
simulations verify our theory. For the case of cubic interactions numerical experiments 
indicate regimes in which only a finite number of species survives, even when the 
thermodynamic limit is considered. 
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1. Introduction 

Replicator equations describe the evolution of self-reproducing interacting species within 
a given framework of limited resources. They have found wide applications in a 
variety of fields including evolutionary game theory, socio-biology, pre-biotic evolution, 
optimization theory and population dynamics, where equivalences to Lotka-Volterra 
equations can be established [H El E]- So-called random replicator models (RRM) 
are also interesting from the point of view of statistical mechanics, as they constitute 
complex disordered systems with non-trivial co-operative behaviour such as phase 
transitions, ergodicity breaking and memory effects. 

Statistical physics can contribute to the analysis of large systems of replicators 
with quenched interactions. Stochastic couplings drawn at random from some fixed 
probability distribution here reflect a lack of knowledge of the detailed interactions of 
real replicator systems. The tools also used for the analysis of disordered systems in 
the physics context are then well suited to compute typical properties of such random 
replicator models, i.e. averages of observables over the distribution of couplings. 
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The first replicator system with quenched random couphngs was introduced and 
studied with path integral methods by Diederich and Opper in lUEHn]. Most subsequent 
studies in the statistical physics community were then based on either replica theory or 
purely on computer simulations [71IHll3linilIIllI21II3IIllISl- 

The restriction to replica methods severely limits the types of models that can 
be addressed, as static analyses require the existence of a global Lypunov (or energy) 
function which is extremised by the replicator dynamics. This in turn requires the 
interaction matrix of the replicators (for pairwise interaction denoted by Jij) to be 
symmetric against permutations of the interacting species, Jij = Jji (generalisation to 
higher-order interactions is straightforward and will be discussed below). This appears 
unrealistic both in the game-theoretic and the biological context. 

In the framework of evolutionary game theory replicator equations describe the 
evolution of one or several populations of players, who engage in a game which is 
played repeatedly The interaction matrices here encode the payoffs given to the 
individual players, these then reproduce according to their success. In the simplest 
case of two players i and j participating in the game Jij would denote the payoff of 
player i when playing against j, and Jji the corresponding payoff for j. Replica theory 
can here be applied to the resulting replicator model only when Jij = Jji for all pairs 
{i,j) of players. This corresponds to equal payoffs for both players, which does not 
appear to be the typical situation in real games. On the contrary one would like to 
introduce an anti-correlation between J^- and Jji so that player i's payoff is positive, 
whenever j's is negative and vice versa. In particular zero-sum games correspond 
to the case Jij = —Jji, i.e. full anti-correlation. Similarly, the {Jij} in biological 
applications describe interactions between different species. The fully symmetric case 
here corresponds to either reciprocal stimulation or reciprocal inhibition, which again 
is unrealistic, as prey-predator relations between two species i and j would for example 
require Jij and Jji to be of opposite signs. 

In order to address replicator systems with fully or partially asymmetric couplings 
methods other than replica theory are required. The only existing analytical studies 
appear to here be the early papers by Rieger |16| and by Opper and Diederich [SJ IH], 
where RRMs were studied using generating functional techniques. 

Such dynamical methods are not concerned with the extremisation of a Lyapunov 
function, but address the temporal evolution of the system directly. Generating 
functionals therefore provide a powerful tool to deal with non-equilibrium disordered 
systems lacking detailed balance. While they were originally developed to study 
magnetic phenomena in spin-glasses [13 CHj, they have recently also been applied to 
models of interacting agents for example in the so-called Minority Game |19| 1^ 1^. 

The aim of the present paper is to discuss dynamical methods in the context of 
random replicator models and to extend the work of Ej (which is concerned with 
pairwise Gaussian interactions) to more general Gaussian replicator systems. To this 
end we work out the effective description of the dynamics of such RRMs in terms of a 
single-species process, and then extract the relevant order parameters in the resulting 
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fixed point states, based on an ergodicity assumption. The breakdown of ergodicity 
in turn signals tfie onset of memory effects, and we compute tlie corresponding pliase 
diagrams. Some differences in tlie types of transitions are liere observed depending on 
tlie symmetry of tlie couplings and the order of the interactions (e.g. pairwise versus 
higher-order). Our study is in direct extension also of [T, where the case of higher-order 
symmetric couplings with Gaussian distribution was addressed using replica theory. 



2. Model definitions 



We consider replicator systems of interacting species labelled by Roman indices 
i = 1,...,N with time-dependent concentrations Xi{t) > 0. Replicator equations are 
generally of the form 

|x.(t)=x.(t)(/.[x(t)]-/[x(t)]), (1) 

with /j[x] the 'fitness' or 'success' of species i (which may depend on the concentrations 
X = (xi, . . . , xn) of all species), /[x] in turn denotes the 'mean' fitness /[x] = J2i ^ifil^]- 
Thus the concentration of species who are fitter than the average increases, and the 
weight of species less fit than average decreases. Note that Eq. ((T)) conserves the total 
concentration J^i^ii'^) — ^■ 

We will here focus on systems with second and higher-order interactions, with 
couplings drawn from a Gaussian distribution, i.e. we will consider replicator equations 
of the form 



^Xi(t) = -Xi{t) 



2uxi{t)+ Jl^,i.,,...,i^Xi^{t)xi^{t)---Xi^{t) - ii{t) -h{t) 



(2) 



with p a fixed integer and where M^^^ = {{i2, . . . ,ip) : 1 < 12 < is < ■ ■ ■ < ip < 
N;i2, . . . ,ip ^ i}. Note the overall negative sign introduced for later convenience. 

The first (diagonal) term in the square brackets describes a self-interaction. Such 
terms were first introduced in the context of replicator models in j2j. Depending on the 
value of the model parameter m > 0, the so-called 'co-operation pressure', the growth 
of any single individual concentration is either strongly or only moderately suppressed 
(large versus small values of u respectively). Given the fixed total concentration of 
species large values of u hence drive the system towards a co-operative state in which 
many species survive at small individual concentrations, whereas small co-operation 
pressure can lead to stationary states in which a large fraction of the species dies out 
asymptotically and where the remaining ones operate at high concentrations. In terms 
of game theory this corresponds to states in which many or few pure strategies are being 
played respectively . In the limit u ^ 00 the only stable attractor is given by a fixed 
point at which all species carry equal weight. 

The second term in the square bracket of (0) describes the interaction between 
species. The couplings {J} are assumed to be drawn from a Gaussian distribution with 
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zero mean and with variance and correlations as follows P^ : 




r 



(3) 



denotes an average over the distribution of the couplings. F G [—1, 1] is a symmetry 
parameter. For F = 1 their distribution is fully symmetric with respect to permutations 
of the indices, F = corresponds to the case of fully uncorrelated couplings, and F = — 1 
to the fully anti-correlated case. Choosing — 1 < F < 1 allows one to interpolate 
smoothly between the different regimes. In the case of pairwise interactions p = 2 the 
choice F = — 1 corresponds to Jij = —Jji, i.e. to zero sum games as mentioned in the 
introduction. For p > 2 the value F = — 1 cannot be realised, since e.g. for p = 3 it 
is impossible to draw three couplings all being equal in their absolute values but with 
pairwise opposite signs. 

The field h{t) in Eq. (0) has been introduced to generate response functions in the 
course of the dynamical analysis of the system. It will be set to zero at the end of the 
calculation. 

Finally, up to a sign the normalisation parameter K{t) corresponds to the mean 
fitness /[x] mentioned above and is chosen to ensure the constraint 



at any time t. We choose initial conditions such that this constraint is fulfilled at the 
starting point t = 0. Note that we have here re-scaled the concentrations such that 
^•Xj = N instead of = 1; this is to ensure appropriate scaling of the statistical 

mechanics quantities and a well-defined thermodynamic limit. 

3. Generating functional analysis and fixed point solutions 

3.1. Effective theory 

Eqs. (121) define a set of N evolution equations, which are coupled to each other through 
the random interactions J, and through the overall constraint (j3)). A Lyapunov function 
governing this dynamics can be found only in the case of symmetric couplings F = 1. In 
the following we will apply generating functionals originally developed by De Dominicis 
[T7] to obtain a macroscopic description of the coupled A^-species dynamics. 
To this end one defines a dynamical generating functional 



((...)) denotes an average over all possible paths of the system. po(x(0)) describes 
the (possibly random) initial conditions from which the dynamics is started. All initial 
concentrations are taken to be strictly positive with probability one here, and we assume 




(4) 
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in the following that the distribution of initial conditions factorizes over individual 
species, i.e. Po(x(0)) = lliPo{xi{0)). 

The general procedure then consists in a computation of the disorder-averaged 
generating functional ^'['0], from which dynamical order parameters can be computed 
as derivatives with respect to the fields {il^i{t), h{t)}. The aim is to formulate a closed set 
of equations describing the temporal evolution of a suitable set of dynamical observables, 
which adequately describe the macroscopic dynamics. The relevant order parameters 
for the RRMs considered here are given by the Lagrange parameter k, = {{n{t))} and 
the correlation and response functions {C, G} in the thermodynamic limit: 

Here (. . .) stands for an average over realisations of the possibly randomly drawn initial 
conditions. Note that t') = for t < t' due to causality. 

The evaluation of the disorder-average in the generating functional can be 
performed along the lines of after a transformation = logXj(t) has been 

performed. The calculation is straightforward, but lengthy. We will therefore not enter 
the details of the mathematical derivation here, but will only report the final outcome. 

One finds that a description of the dynamics can be obtained in terms of an effective 
process of the form 



^x{t) = -x{t) 



2ux{t)-r^ — y / dt'G{t,t')C{t,t'y~^x{t')+r]{t)-K{t)-h{t) .(7) 

2 .In 



This equation describes a stochastic process for a representative species concentration 
{x{t)}, which is a random variable and subject to Gaussian noise {f]{t)}, with non-trivial 
temporal correlations according to 

mvit')). = lc{t,t'r'. (8) 

The correlation and response functions C and G of the original multi-species system 
(see Eq. (jH}) can in turn be shown to be given by 

C{t,t') = {x{t)x{t')l, G{t,t') = (^^'^ , (9) 

and K in (jZj) has to be chosen such that 

(x(t)), = 1 Vt. (10) 

(. . .)^ denotes an average over reahsations of the effective process (jlj), i.e. over 
reahsations of the single-species coloured noise {ri{t)} and over initial conditions 
described by po{x{0)). 

Eqs. ()7l8l9ll0p thus form an implicit, but closed set of equations from which 
{k, C, G} are to be obtained. For p = 2 this system coincides with the results of 0. 
The description in terms of an ensemble of decoupled, but stochastic effective species 
is exact and fully equivalent to the original iV-species problem (defined by Eq. (0)) 
in the limit — > cxo (in the sense that disorder-averages of macroscopic observables 
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in the original dynamics can be obtained as averages (• • •)^ on the level of the effective 
process) . The retarded self-interaction and the coloured noise in the effective process are 
direct consequences of the quenched disorder in the original problem and impede a full 
explicit analytical solution of the self-consistent saddle-point problem for the two-time 
quantities C{t,t') and G{t,t') and the full function K{t). One is therefore limited to a 
specific ansatz for the trajectories of the effective particles. 



3.2. Fixed point solutions 

In order to proceed analytically we assume that the system reaches a fixed point 
asymptotically, i.e. Xi{t) — Xj as t — >^ oo for all i in the original dynamics. Numerical 
simulations confirm that this assumption is indeed justified for values of u larger than 
some critical value Uc (the value of Uc may depend on p and T as discussed in detail 
below) . 

We may then make a similar assumption for the realisations of the effective process, 
i.e. x{t) — i> X, with x a static random variable. Accordingly the correlation function 
becomes flat at the fixed point and K,{t) approaches a fixed point as well, so that we 
write 

lim C{t + T,t) = q Vr, lim K{t) = k. (11) 

>oo t— >oo 

We also make the standard ansatz of a time-translation invariant asymptotic response 
function 

lim G(t + r,t) = G'(r). (12) 

t— >oo 

Furthermore we will only address ergodic stationary states, that is states in which 
perturbations have no long-term effects so that the integrated response function remains 
finite 

I drGir) <oo, (13) 
Jo 

and no long-term memory is present, i.e. we will assume 

limG(t,t')=0 (14) 

also for finite t'. In the absence of memory we may send the starting point of the 
dynamics to — cxd for convenience. 

Finally within the fixed-point ansatz we also assume that each realisation of the 
the single-species noise {rjit)} approaches a time-independent value rj asymptotically, 
which according to (jH)) is then a static Gaussian variable with zero mean and variance 

(^^>. = (15) 
Fixed points of ((7j) then fulfill the condition 

X Uux - V ^^^~^\ p-\x + (|/"') ^'^ z-t^=Q, (16) 



X 
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1 /2 

where we have written rj = z with z a standard Gaussian variable. h{t) 

has been set to zero at this stage (derivatives with respect to h{t) can be obtained 
as derivatives with respect to the effective-particle noise r]{t) up to an inverted sign). 
We note that x{z) = is a solution of this equation for all realisations of the random 
variable z. Other solutions may be possible whenever setting the bracket to zero leads 
to a positive value of x{z). Taking into account the stability of the zero solutions, one 
finds that the ansatz 



r z) = — - B 

2u-rp{p-l)xqP-y2 



2^ ) ' 



(17) 



adequately describes the solutions, with 9[x] the step function (9[x] = 1 for x > 
and 9[x] = otherwise). Similarly to EH] one demonstrates that zero fixed points 
x{z) = are unstable for positive arguments of the 9-function. Self-consistency then 
demands that 

Dzx{z) = l, [Dzx{zf = q, - (^q^'-^Y'^^ f Dz^4^ = x, (18) 



dz 

with Dz = ^e-" /2 ^ 

standard Gaussian measure. Upon using the explicit ansatz p7|l 
for the fixed points these equations may be written as 



A 



^q'-Y^'(2u-rP^P^Xq'-') = j ^Dz{A-z), (19) 
q ( ^q^-') 'U2u- rP^^xq'-'] ' = ^ Dz{A - z)\ (20) 



2 



oo 
A 



2u - r^^^x/-^ ) X = / Dz, (21) 



with A = k/ {^q^'^Y^"^ ■ We note that (f) = j^^Dz represents the probability of a 
given species i to survive in the long-term limit, i.e. to attain a fixed point value 
Xi = limt^oo Xi{t) > 0. We will refer to as the fraction of surviving species in the 
following. For the case of symmetric couplings, F = 1, equations EDI EH) are 
identical to those found in replica-symmetric studies of the statics of the model [7], for 
p = 2 they furthermore coincide with the ones reported in [Sj. They are readily solved 
numerically (in terms of k, q and x) for arbitrary values of the model parameters u 
and F (at fixed p)|. Results for the order parameter q are depicted for the model with 
pairwise interaction {p = 2) in Fig. ^ We find near perfect agreement between the 
analytical theory and numerical simulations§ for u larger than some critical value Mc(F) 
(for F < Tc{u) respectively in the right panel of Fig. QJ. The deviations from the theory 
below Uc (resp. above Fc) are due to a breakdown of some of our assumptions regarding 

t Note here that the right-hand-sides of dl EOI EB) 

can be written in closed form as functions of A 
after performing the Gaussian integrals over z. A solution of these equations can thus most efficiently 
be obtained by expressing {u, q, x} as functions of A, and by subsequently varying A 
§ A discretisation method similar to that of [H] is used. The typical (effective) time-step is of the order 
of Ai« 0.01 -0.1. 
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Figure 1. Inverse order parameter q for the model with p = 2. Symbols are 
from simulations for N — 300 species, run for 20000 discretisation steps (effective 
time-stepping is of order 0.01 — 0.1) and averaged over 50 samples of the disorder. 
Measurements are taken in the last quarter of the simulations. Initial conditions 
correspond to a flat distribution over the interval [0,2]. The solid lines represent 
the analytical predictions in the ergodic phase, and have been continued as dashed 
lines into the region of broken ergodicity, where our theory is no longer valid. Left: 
1/q as a function of u for F = —1,-0.5,0,0.5,1 from top to bottom. Right: Same 
quantity as a function of F at fixed w = 2, 1, 0.5 (top to bottom). 



ergodicity and the stability of the fixed points, the resulting phase transitions will be 
discussed in more detail below. 

Let us finally in this section turn to an interpretation of the order parameter q, 
which within the fixed point ansatz is given by g = lim^r^oo (^^f)- If the Xi 

were normalised to 1 instead of N, q would be closely related to what is known as 
Simpson's index of diversity [22] and would indicate the probability that two randomly 
selected individuals of the eco-system are of the same species. A Simpson index of 
zero thus indicates infinite diversity, a value of one corresponds to no diversity. A 
similar interpretation holds for the normalisation of the {xi} used here: finite values 
of q indicate the co-existence of an extensive number of species with concentrations 
Xi ~ 0{N^), while a divergence in q signals the dominance of a sub-extensive number 
of species with diverging concentrations in the thermodynamic limit jH|. 1/g is thus a 
measure for the diversity of the system, and we expect the behaviour of 1/g to be similar 
to that of the fraction of surviving species = / Dz Q[x{z)]. This is indeed confirmed 
in numerical simulations (not shown here) in which one verifies that 1/q and show 
the same qualitative behaviour as functions of the model parameters. The left panel of 
Fig. n]thus confirms the role of m as a co-operation pressure, the higher the value of u 
the more species survive in the long-term limit and the more diverse the system is in 
its stationary state. The infiuence of the symmetry of the interactions on the diversity 
(right panel of Fig. [TJ will be discussed below. 
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3.3. Stability analysis and phase transitions 

Our theory and solution crucially relies on the assumptions made with respect to 
ergodicity, weak long-term memory, and the stability of the fixed point which the 
dynamics is assumed to approach asymptotically. A breakdown of this theory may 
thus be signalled by the failure of any of these assumptions, marking a dynamical phase 
transition at which the ergodic fixed-point regime ceases to exist. In general one may 
imagine the following phenomena to occur: 

(i) The analytical ergodic theory resulting in Eqs. (|T8|) might predict the integrated 
response x to diverge in some subset of the parameter space of the model. Such a 
transition has been observed in replicator models with Hebbian interactions j28j . 
and below this transition the system is found to be non-ergodic. 

(ii) In a similar way the theory might analytically predict singularities in either q or k. 
No such types of transitions seem to have been observed so far in studies of RRM||. 

(iii) The set of equations ()18j) might fail to have solutions for certain values of the 
control parameters, so that no ergodic states are allowed in this region of the phase 
diagram. We will discuss cases of this type in section HJ 

(iv) Our ergodic theory is based on the assumption that all trajectories of the system 
will evolve into fixed points asymptotically. An onset of instability of these fixed 
points against small perturbations will thus signal the breakdown of our theory as 
they will no longer be local attractors of the dynamics. Such transitions have been 
observed in RRM with Gaussian and Hebbian couplings in El l2Hj • 

(v) Finally a breakdown of our requirement that long term memory be weak would 
indicate that the system remembers perturbations during the transient dynamics. 
Even if the system still evolves into a fixed point the latter might no longer be 
unique and the choice of the asymptotic stationary state might thus depend on 
initial conditions or perturbations at early stages of the temporal evolution. This 
type of transition has been related to the previous one (instability of the assumed 
fixed point) in |2HI; and we will demonstrate below that the conditions for both types 
of transitions are indeed fulfilled at the same points in parameter space whenever 
a replicator system with symmetric couplings is considered. 

Transitions of the types (i)-(iii) are easily detected in numerical solutions of Eqs. 
fll9l20l21|) (or their analogues for other types of RRMs) . In the following we will therefore 
focus on (iv) and (v) and will derive explicit conditions for the onsets of instability and 
of memory at finite integrated response respectively. 

3.3.1. Instability of the fixed point In order to inspect the stability of the assumed 
fixed point we will follow 5J and add a small component sC,{t) of white noise to the 

II Phases in which hiiiAr^oo = have been hinted at in and will be discussed in more detail 
below. Note however that in these cases the equations describing the ergodic states do not analytically 
lead to singularities in but do not allow for any solutions at all in the phases where q is found to 
diverge in simulations; they are of the type (iii) in the above list of possible transitions. 
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system (with {({t)) = 0, {((t)((t')) = 6(t — t')) and then study small fluctuations ey(t) 
and ev{t) of the effective species concentration and the single-particle noise about their 
respective flxed points x and rj: 

x{t) = x + ey{t), 7]{t) = r] + ev{t). (22) 

In the following we will only consider the case x > 0, as it turns out that the onset 
of instability in replicator systems is determined by the flxed points different from 
zero while the zero flxed points remain stable against perturbations (see also [3| I28j). 
Inserting (j22I) into the effective process ((Zj) (for flxed points x > 0) and adding the 
additional noise component leads to 



d , , 



2uy{t) - r^^^-ilg^'-^ f dt'G{t - t')y{t') + v{t) + C(t) 



(23) 



to flrst order. In Fourier space one has 



yioj) = — : — -T . (24) 

f + 2m - r2^gp-2G(cu) 

(with {y{u),v{u!)X{^),G{u!)} the Fourier transforms of the one-time functions 
{y{t), v{t), C(t), G(t)}). Focusing on u = one obtains 

(|y(^ = O)P>^ = 0x Kl^(^ = 0)O. + l]^. (25) 



The flrst factor of in (j25j) takes into account that ()24jl was derived for non-zero 
flxed points a; > and that fluctuations about zero flxed points do not contribute to 
{\y{uj = 0)p)^. We have also used Eq. (ED) to replace (2m - Tp{p - l)qP-\/2)-'^ by 
/ cfp'. The self-consistency condition on the covariance of the single-particle noise (Eq. 
(jH))) on the other hand implies 

(R-)r>. = ^^/^^(iyMp>. (26) 

Insertion of (pUj) into then allows one to solve for (|y(0)p)^, and one flnds 

.X 2 

We conclude that (|y(0)p)^ diverges when the square bracket on the right-hand-side 
of (jTfjl vanishes, suggesting an onset of instability. (|5'(0)P)^ is predicted to become 
negative whenever the right-hand-side of (P7j) becomes negative, indicating a further 
contradiction. The onset of instability thus occurs at the point deflned by 

'-^^^^X^ = 1. (28) 

and the flxed points are stable whenever this expression is strictly smaller than one. 
With some further algebra one can show that ()28|) deflnes a line UciT^p) = Uc{T = 
0,p)(l + r) in the {u, F) plane for any flxed value of p. 
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3.3.2. Memory onset In the previous section we have related the breakdown of our 
ergodic theory to a local instability of the fixed point reached by the dynamics. We 
will now inspect for an onset of long-term memory at finite integrated response. This 
type of transition has been observed previously for example in Minority Games with 
self- impact correction or diluted interactions |2S1 and can also be interpreted in 
terms of a breakdown of time-translation invariance, see [201 for details. It is also found 
in replicator systems with Hebbian interactions [2H]- In order to see how solutions with 
memory bifurcate from the time-translation invariant ergodic states we will proceed 
along the lines of [201 and make the following ansatz for the response function 

G{t,t') =Go{t-t')+eG{t,t'), (29) 

where £G(t, t') is a small contribution which breaks time-translation invariance. The 
starting point of the dynamics can here no longer be sent to —oo. For a physical 
interpretation of Go and G see |2H]- Following [20] G is taken to depend only on the 
(earlier) time t' asymptotically, i.e. lim(_,oo G'(t, if:') = G{t'). 

Starting from one expands the kernel of the retarded self-interaction in the 
effective process to linear order in G, and finds the following effective process: 

^xit) = - xit) huxit) - r^fc^ f dt'Goit - t')C(t - ty-^xif) (30) 

at L 2 Jo 

_^p P(p-l) I dt'G{t')C{t-t'y-^x{t') + 7]{t)-K{t)], (31) 
so that non-zero asymptotic fixed points are now found to fulfill 

2ux - T ^^^~^\ qP-'x - eV ^^P'^^ f dt'G{t')q^~^x{t') + C^q''^) "\-k = Q. (32) 
2 2 Jq V2 / 

From this one can compute the response function G{t") at a transient time t" self- 
consistently and finds 

Git") = TP^f-^(2u-TP-^xf~']" fdt'Git')/^) . (33) 



V - y JO 

This can be understood as an eig envalue problem of the form G{t) = J dt' M{t-t')G{t'), 
with a suitable kernel M [20] • Upon taking Fourier transforms and focusing on the zero 
frequency mode = one finds after some simpifications 

X = r^^^^-^x^x, (34) 

with X = / dtG{t). Although x = is always a solution, one realises that non-zero 
solutions for £ become possible at the point at which 

^Pip^^. = 1 (35) 
20 ^ ^ ' 

We will refer to this as the memory onset (MO) condition in the following. 
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Figure 2. Phase diagrams for the models wit p = 2 and p = 3 in the (tt, F) plane. Left 
panel {p — 2): the dashed line marks the onset of long-term memory, see Eq. Ij^SII . the 
solid line marks the onset of the instability of the fixed points, Eq. (|28|l . Right panel 
(p = 3): the solid curve separates the region in which Eqs. H19|l - H21() admit a solution 
(high u) from the region in which no solutions are found, dashed line marks the onset 
of instability, Eq. ISHJ- 



4. Phase diagram and discussion of results 

We now turn to a discussion of the resulting phase diagrams, as depicted in Fig. |21 In 
general we find that our ergodic fixed-point ansatz is valid and stable for values u > Uc 
with Uc some critical value which depends on the order of interaction p and on the 
symmetry parameter F. Above Uc the system evolves into a unique fixed point, with 
virtually all species surviving in the limit u ^ oo. This corresponds to a fixed point in 
the interior of the simplex defined by the condition A^~^ Yli^i ~ ^- -^^^ small values of 
u fewer species survive so that the system operates near the boundary of this simplex. 

The dependence of the phase diagram on F is particularly interesting. One notes 
that both conditions (PHj) and (jH^I) coincide for symmetric couplings F = 1, but that 
the onset of the instability of the fixed point occurs before memory sets in for F < 1 ^. 
At the same time simulations reveal different types of behaviour below the transition, 
depending on the symmetry of the interactions: while for fully symmetric interactions 
F = 1 the system reaches a fixed point also below Uc , fixed points are not necessarily 
observed below Uc for even partially asymmetric couplings. This will be discussed in 
more detail below. The type of transition also depends on whether the interactions are 
pairwise or cubic, we find that in the system with p = 2 the number of surviving species 
is always extensive for any u > 0, whereas this number can turn out to be 0{1) in 
the model with p = 3 below the transition even in the thermodjTiamic limit. In the 
following we will treat the cases p = 2 and p = 3 separately. 

% The equivalence of the MO and instability conditions at F = 1 can be extended to more general 
replicator equations derived from a Lyapunov function. 




Figure 3. Relative distance <P /q and roughness h/q versus u for the model with 
p — 2. Connected symbols are from simulations (triangle up: V = —1, triangle left: 
r — —0.5, circles: F = 0, squares: F = 0.5, diamonds: F = 1). is obtained by 
comparing the stationary states of two runs at fixed disorder started from independent 
random initial conditions over [0,2]. Simulation parameters as in Fig. n Vertical 
dashed lines indicate the location of the phase transitions as predicted by the theory 
for F = — 0.5,0,0.5,1 from left to right. No transition occurs for w > and F = — 1. 



4-1. Pairwise interaction (p = 2) 

We find here that Eqs. 1)191201211) admit solutions for all u and F and we observe no 
singularities in any of the persistent order parameters. The predictions of the analytical 
theory for the order parameters in the ergodic phase are compared with numerical 
simulations in Fig. ^and we find near perfect agreement in the ergodic phases {u > Mc(F) 
and F < Fc(m) respectively). 

The only possible violations of our assumptions of the ergodic stationary state are 
instabilities of the fixed points or an onset of long-term memory at finite susceptibility 
X- The condition for onset of the instability of the fixed points fl28|) reduces to 
Uc{p = 2, F) = (1 + F)/(2v^), as already derived in [5j. If one ignores this instability 
and continues the solutions of ()19I2UI2H) to the right of the line u = Uc{p = 2, F) one 
finds that the memory onset condition Eq. ()35|) is fulfilled along the dashed line in the 
left panel of Fig. |21 As the fixed points are already unstable in this regime, this line 
has no direct physical meaning though, and is depicted only for illustration. We note 
that for F = 1 both the instability of the fixed point and the onset of memory occur at 
the same value of m = l/V^. In the statics a de Almeida-Thouless instability is found 
at this point and replica symmetry is broken at smaller values of u jT^. Interestingly, 
one does not find any transition at any positive co-operation pressure u > for full 
anti-correlation F = —1, i.e. in the case of zero-sum games. 

To verify the onset of non-ergodicity further we have performed simulations in which 
two copies of the system with the same realization of the disorder are generated and in 
which the dynamics is started from two independent sets of random initial conditions 




Figure 4. Trajectories {xi{t)}t for 10 individual species in single runs of simulations 
of the model with p = 2 at fully asymmetric couplings, F = 0. Left: u — 0.2 below 
the transition, where fixed points are unstable; right: u — 1.0 above the transition. 
Simulations are for N — 500 agents runs for 20000 time steps. 



drawn from a flat distribution over Xj(0) G [0,2]. The two resulting trajectories of 
the system are labelled by {x^it)} and {x'lit)} respectively and we have measured the 
squared distance (f = N~'^ ~ (^))^ the stationary state. Results are 

depicted in Fig. El (left panel). For convenience we plot the re-scaled distance (f/q. We 
find that d"^ is essentially zero above Uc{T), confirming the ergodicity and independence 
of the stationary state on initial conditions in this regime. Below Uc we find that d"^ > 0, 
indicating the existence of multiple stationary states and the pronounced relevance of 
initial conditions. For F = — 1 no memory effects are observed at any m > as expected 
from the theory. 

We have secondly measured the relative 'roughness' h/q of the trajectories, where 
h = J^iii^iity) t~ {^ii^))^) (• • denotes a time-average in the stationary state. 
As shown in the right panel of Fig. IHl all trajectories are fiat [h = 0) irrespective of F 
above Uc(F), verifying that the system indeed evolves into a fixed point in this regime. 
For symmetric couplings F = 1 fixed points are also reached below the transition. For 
— 1 < F < 1 we find non-zero values of h below the transition, so that we conclude that 
the system does not necessarily evolve into fixed points here"'". This behaviour is also 
confirmed in Fig. E] where we show typical trajectories of the system at F = (fully 
uncorrelated couplings) below and above and the transition, see also jHI for a similar 
figure. For F = 1 one finds that the system evolves into a fixed point irrespective of u 
(trajectories not shown here). 

We conclude the section concerning the model with p = 2 by some brief summary of 

+ The numerical data for h appears to show some dependence on system size, running time and time- 
window over which measurements are taken in the regime of volatile trajectories below Uc{T < 1) (i.e. 
when h > 0). The observations that fixed points are attained for all F above Uc and for all u for F = 1 
are however not affected. More extensive simulations seem appropriate in order to check whether the 
apparently decreasing values of /i as u approaches zero at F < 1 is genuine or a numerical artifact. 
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Figure 5. Average fitness / = —k versus F for tlie model with p = 2. Symbols 
are from simulations {u = 0.5, 1,2 from top to bottom), simulation parameters as in 
Fig. The solid lines marks the predictions of the analytical theory, for u = 0.5 the 
theoretical line has been continued as dashed line into the phase where the ergodic 
theory is no longer valid. 

the effects of asymmetry in tlie couplings. As discussed above even partial asymmetry 
changes the nature of the phase transition as compared to the fully symmetric case. This 
was already pointed out in [SJ . For symmetric couplings the system evolves into a fixed 
point asymptotically both above and below Uc(T = 1) = 1/V2. Below the transition 
fixed points are however not unique, hence the breaking of ergodicity and dependence on 
initial conditions. For F < 1 we find that the system reaches a fixed point asymptotically 
above Mc(F), but not necessarily below, where instead one finds volatile trajectories in 
which species can almost die out and then recover to large concentrations at later times. 
No transition is observed for F = — 1 (full anti-correlation), and the system here evolves 
into a unique fixed point for all m > 0. The effects of asymmetry on the diversity 
1/q are demonstrated in Fig. ^ (right panel). Generally speaking asymmetry increases 
the diversity of species in the stationary state (equivalently decreasing F leads to a 
higher fraction of surviving species) . While this effect is small for large values of u some 
significant dependence on F is observed for small co-operation pressure u. 

Finally the effects of asymmetry on the fitness of the species in the stationary state 
are demonstrated in Fig. El The average fitness is here given by / = Yli^ifii^] 
with /j[x] = —2uXi — JijXj, and it turns out that fi = f for all surviving species 
i. From the analytical calculation one finds / = —k. As shown in the figure an 
increased asymmetry in the couplings leads to a reduced overall fitness. Again the 
effect is strong for small u but negligible for large co-operation pressure (where the 
diagonal term proportional to u dominates the interactions). Thus, in systems with 
symmetric couplings fewer species tend to survive than in systems with uncorrelated 
or anti-correlated interactions, but these surviving species have higher fitness than the 
many survivors of systems with asymmetric or anti-symmetric couplings. 
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Figure 6. Reciprocal order parameter q and distance (P /q versus u for the model 
withp = 3. Symbols are from simulations (triangles: F = —0.5, circles: F = 0, squares: 
F = 0.5, diamonds: F = 1). Solid lines for 1/q are the analytical predictions in the 
ergodic phase u > Uc, continued as dashed lines into the non-ergodic phases. Vertical 
lines indicate the predicted locations of the transitions. Simulations were performed for 
N = 200 species, started from random initial conditions drawn with flat distribution 
from Xi{0) e [0,2] and run for 3000 time-steps. Results are averages over 20 samples 
of the disorder. 



4-2. Cubic interaction (p = 3) 

For p = 3 we find a distinctively different behaviour of tlie system. Eqs. fjl9l QUI I7H) 
admit solutions only for u > u*{r), and below this value no ergodic fixed point solutions 
can be found. The obtained value u*{r = 1) ^ 2.028 has already been reported in the 
context of a replica analysis of the model with symmetric couplings [Tj. Note also that 
the limits lim^jj^j. q, lim„j„* x and lim^j^u* k exist and take finite values. As indicated 
in Fig. 121 this type of transition is preceded by an onset of instability of the fixed 
points (determined by Eq. (j28|) ) for values of F < Fq ~ 0.7 in the sense that the 
instability sets in as u is lowered starting in the ergodic phase before solutions of the 
saddle point equations cease to exist. We note that the onset of instability again occurs 
along a line in the (m, F) plane, with Uc{p = 3, F) = Uc{p = 3, F = 0)(1 + F), where 
Uc{p = 3, F = 0) ~ 1.028*. We find that solutions of the saddle-point equations exist 
for all u > if F is (sufficiently) negative jj- 

* It turns out that Eqs. H19I20I21|) have two branches of solutions above u*, both solutions merge at 
u* . For F > Fo we find that the fixed points are stable along the physical branch and that the condition 
(I28II for the onset of instability is fulfilled at Uc{p = 3,F) — Uc{p = 3,F = 0)(1 + F) only along the 
unphysical branch. Thus no instability is observed and a transition of the above type (iii) is the first 
one to occur when u is lowered from infinity. For F < Fq the fixed point becomes unstable along the 
physical branch and the instability transitions precedes the failure of solutions to exist, 
tt Some numerical inaccuracies are encountered for F « when determining the line along which 
solutions cease to exist. The solid line in the phase diagram for p = 3 (right panel of Fig. could 
therefore be investigated more carefully near F = 0. At this point instabilities of the fixed points have 
already set in and solutions are unphysical, so that the precise location of the point at which solutions 
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Figure 7. Left: Exponents 7 (circles) and ~6 (diamonds) (see text) as functions 
of u for the model with p — 3 {T = 1). Each data point is obtained by performing 
simulations at fixed u for = 50, 75, ... , 200, simulations are run for 6000 steps, 
averages over 10 samples of the disorder are taken, and subsequent fits to power laws 
q ~ and (p ~ are performed to obtain 7 and S. Vertical bars indicate the 
standard error of 7. Vertical solid line marks u*{T — 1). Inset: Log- log plots oi q vs N 
for u = 1.0 (circles), u = 1.5 (squares) and u = 2.5 (diamonds) and the corresponding 
fits to power laws. (The data for q at u = 2.5 in the inset has been multiplied by a factor 
of two for convenience.) Right: Temporal evolution of the number iV+ of species with 
concentration Xi{t) > 0.01 for the same parameters, and at u = 1. Different curves are 
obtained from simulations of different system sizes N = 50, 100, 150, 200 from bottom 
to top. 



Results from numerical simulations of the model with p = 3 are presented in 
Fig. |B1 Given the cubic interaction we limit ourselves to relatively small system sizes 
= 200 and to moderate running times here. Results in the ergodic phases are however 
not affected significantly by these constraints in computing time. We observe good 
agreement between the results for the order parameter q as obtained from theory and 
experiment in the ergodic region, with only slight discrepancies close to the transitions 
due to finite-size effects. Again asymmetry in the couplings increases the diversity and 
number of surviving species. The right panel of Fig. IHl confirms that non-ergodicity 
sets in at Uc{T), indicated by non-zero distances between the stationary states obtained 
when starting the dynamics at fixed disorder from different random initial conditions as 
explained above. 

To summarise we have found a regime u > max{uc{T),u*(r)) in which our ergodic 
theory applies and matches the simulations. For F < Fq we then find an intermediate 
regime < m < in which solutions of the saddle point equations exist, but in 
which the attained fixed points are unstable and the stationary state depends on initial 
conditions. No such intermediate regime appears to be found for F > Fq. Finally below 
u* no solutions of the equations for the order parameters in the stationary ergodic state 
exist. The transition at u* is not marked by any divergence in the analytic solutions 
for {q,x,i^) as u I u*. To conclude this section we report some indications about 

cease to exist is immaterial for the behaviour of the model. 
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the behaviour of the system below u* based on numerical simulations (we restrict to 
r = 1). While the persistent order parameters converge to finite values in the limit 
of large above u* we observe strong finite size effects for u < u* . The numerical 
data are consistent with power law behaviour q ~ A^'*' and ~ A^"^ with 7 > and 
5 < below u* , see Fig. |3 This divergence of q with seems to suggest that only 
a sub-extensive number of species (j)N ~ A^^+'' survives in this regime [5 < 0) in the 
thermodynamic limit. Given the constraint limjy^oo Tlii^i ~ ^1 then conclude 
that these surviving species each must have concentrations scaling as Xi ~ A^""^, which in 
turn results in g = A^"-*^ ~ A^"'', i.e. we expect 7 = —5. Within the accuracy of our 

simulations this is indeed what we observe, see Fig. More precisely our simulations 
are consistent with 7 = —5 = 1 at low values of m, so that we expect only a finite number 
of species to survive in this regime {(pN ~ 0{N^)), even in the thermodynamic limit. 
This unusual effect is indeed confirmed upon plotting the number of species which are 
still alive as a function of time for different system sizes, see the right panel of Fig. [7| 
Closer inspection shows that the final number of surviving species in these simulations 
all fall into the range 10 < Ncj) < 12, even if the system size is varied by a factor of four 
from A^ = 50 to A^ = 200. A similar effect was observed in a replicator model with an 
explicit extinction threshold ^3]. 

While our simulations confirm convincingly that 7 = —5 = above u*, and while 
they seem to suggest that 7 = —6 = 1 sufficiently far below the transition, the accuracy 
of our numerical experiments do not allow us to make any statement as to whether 
this onset of non-extensivity occurs continuously or discontinuously, i.e. whether the 
exponents drop instantly to zero when the transition at u = u* is crossed (from below) 
or whether one observes a cross-over. 

5. Concluding remarks 

In summary we have demonstrated how generating functionals can be used to study 
random replicator models. In the case of symmetric couplings this approach is able 
to reproduce the equations describing the ergodic stationary states obtained from 
static replica studies, but in addition the analysis of the dynamics also allows to 
address replicator systems with asymmetric interactions and to make accurate analytical 
predictions for the order parameters in their fixed point regimes and the resulting phase 
diagrams. 

Our analysis of Gaussian RRM indicates that the effect of asymmetry and anti- 
symmetry in the couplings is to increase the stationary value of the diversity parameter 
1/q and with it the number of surviving species. This is the case for both models studied 
here, with interactions between p = 2 and p = 3 species respectively. Asymmetry in the 
couplings also reduces the fitness of the surviving species. Furthermore the range of the 
ergodic phases seem to be consistently larger the higher the degree of anti-symmetry in 
the couplings becomes. The symmetry or otherwise of the interactions also determines 
the behaviour of the system in the non-ergodic phases. Replicators with symmetric 



Random replicators with asymmetric couplings 



19 



couplings evolve into fixed points in all regions of the phase diagram. In the symmetric 
case one expects a large number of such fixed points in the phase below Uc and 
initial conditions determine which of these is reached, as indicated by the memory- 
onset. No such phase (and hence no MO transition) is found in models with asymmetric 
or anti-symmetric couplings. Instead, volatile trajectories are here observed below the 
transition, and individual species may come close to extinction Xi{t) -C 1, but then 
recover to macroscopic values Xi{t) ~ 0{1). 

Models in which the self-interaction is of a lower order than the terms coupling 
distinct species appear to exhibit a phase in which the fraction of surviving species 
scales as ~ so that only a finite number of species survives in the long-run 

even in the thermodynamic limit. Possibly parallels with condensation effects in growth 
models of complex networks with quenched fitnesses of links and/or nodes can be drawn, 
as the dynamics of such networks obeys growth laws which are similar to the replicator 
equations studied here ISO]. 

While we have here concentrated on relatively simple single-population systems 
with Gaussian couplings, similar approaches may be taken to address more intricate 
models of replicators. These could include multi-population models of game theory 
or dynamical models of chemical metabolic reactions, which progress in proportion to 
the concentration of reaction partners and in which stoichiometric coefficients might be 
modelled in terms of quenched random couplings. Multi-population replicator equations 
of game theory correspond to so-called bi- matrix games [2], and it would here be 
interesting to study the relation between the fixed points of the replicator dynamics 
and the Nash equilibria of such games jHI], also as a function of the degree of anti- 
correlation in the payoff matrices and the co-operation pressure. Some work is currently 
in progress along these lines. 
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